*DECK CQCBY
      SUBROUTINE CQCBY (LUN, KPRINT, IPASS)
C***BEGIN PROLOGUE  CQCBY
C***SUBSIDIARY
C***PURPOSE  Quick check for SLATEC subroutine
C            CBESY
C***LIBRARY   SLATEC
C***CATEGORY  C10A4
C***TYPE      COMPLEX (CQCBY-C, ZQCBY-Z)
C***KEYWORDS  QUICK CHECK, CBESY
C***AUTHOR  Amos, Don, (SNL)
C           Goudy, Sue, (SNL)
C           Walton, Lee, (SNL)
C***DESCRIPTION
C
C *Usage:
C
C        INTEGER  LUN, KPRINT, IPASS
C
C        CALL CQCBY (LUN, KPRINT, IPASS)
C
C *Arguments:
C
C     LUN    :IN  is the unit number to which output is to be written.
C
C     KPRINT :IN  controls the amount of output, as specified in the
C                 SLATEC Guidelines.
C
C     IPASS  :OUT indicates whether the test passed or failed.
C                 A value of one is good, indicating no failures.
C
C *Description:
C
C   CQCBY is a quick check routine for the complex Y Bessel function
C    generated by subroutine CBESY.
C
C   CQCBY generates sequences of Y Bessel functions from CBESY
C    and checks them against the evaluation from the formula
C
C    Y(FNU,Z*ROT) = C(FNU+1)*I(FNU,Z)-(2/PI)*CONJG(C(FNU))*K(FNU,Z)
C
C    where ROT = EXP(PI*I/2) , C(FNU)=EXP(PI*FNU*I/2) , I**2=-1
C
C    and -PI.LT.ARG(Z).LE.PI/2, in the (Z,FNU) space.
C
C***REFERENCES  Abramowitz, M. and Stegun, I. A., Handbook
C                 of Mathematical Functions, Dover Publications,
C                 New York, 1964.
C               Amos, D. E., A Subroutine Package for Bessel
C                 Functions of a Complex Argument and Nonnegative
C                 Order, SAND85-1018, May, 1985.
C***ROUTINES CALLED  CBESI, CBESK, CBESY, I1MACH, R1MACH
C***REVISION HISTORY  (YYMMDD)
C   830501  DATE WRITTEN
C   890831  Revised to meet new SLATEC standards
C
C***END PROLOGUE  CQCBY
C
C*Internal Notes:
C   Machine constants are defined by functions I1MACH and R1MACH.
C
C   The parameter MQC can have values 1 (the default) for a faster,
C   less definitive test or 2 for a slower, more definitive test.
C
C**End
C
C  Set test complexity parameter.
C
      INTEGER  MQC
      PARAMETER (MQC=1)
C
C  Declare arguments.
C
      INTEGER  LUN, KPRINT, IPASS
C
C  Declare external functions.
C
      INTEGER  I1MACH
      REAL  R1MACH
      EXTERNAL  I1MACH, R1MACH
C
C  Declare local variables.
C
      COMPLEX  CI, CIP, COE1, COE2, CSGN, CSPN, CW, CWRK, V, W, Y, Z, ZN
      REAL  AA, AB, AER, ALIM, ARG, ATOL, AV, CC, CT, DIG, ELIM, EPS,
     *   ER, ERTOL, FFNU, FILM, FNU, FNUL, HPI, PI, R, RHPI, RL, RM,
     *   R1M4, R1M5, R2, SLAK, ST, T, TOL, TS, XN, XNU, XX, YN, YY
      INTEGER  I, ICASE, IERR, IFNU, IL, IR, IRB, IT, ITL, I4, K, KDO,
     *   KEPS, KK, KODE, K1, K2, LFLG, MFLG, N, NL, NU, NUL, NZ, NZ1,
     *   NZ2
      DIMENSION  AER(20), CIP(4), CWRK(20), KDO(20), KEPS(20), T(20),
     *   V(20), W(20), XNU(20), Y(20)
      DATA CIP(1), CIP(2), CIP(3), CIP(4) /
     *   (1.0E0,0.0E0), (0.0E0,1.0E0), (-1.0E0,0.0E0), (0.0E0,-1.0E0) /
C
C***FIRST EXECUTABLE STATEMENT  CQCBY
      IF (KPRINT.GE.2) THEN
        WRITE (LUN,99999)
99999   FORMAT (' QUICK CHECK ROUTINE FOR THE Y BESSEL FUNCTION FROM ',
     *    'CBESY'/)
      ENDIF
C-----------------------------------------------------------------------
C     Set parameters related to machine constants.
C     TOL is the approximate unit roundoff limited to 1.0E-18.
C     ELIM is the approximate exponential over- and underflow limit.
C     exp(-ELIM).lt.exp(-ALIM)=exp(-ELIM)/TOL    and
C     exp(ELIM).gt.exp(ALIM)=exp(ELIM)*TOL       are intervals near
C     underflow and overflow limits where scaled arithmetic is done.
C     RL is the lower boundary of the asymptotic expansion for large Z.
C     DIG = number of base 10 digits in TOL = 10**(-DIG).
C     FNUL is the lower boundary of the asymptotic series for large FNU.
C-----------------------------------------------------------------------
      R1M4 = R1MACH(4)
      TOL = MAX(R1M4,1.0E-18)
      ATOL = 100.0E0*TOL
      AA = -LOG10(R1M4)
      K1 = I1MACH(12)
      K2 = I1MACH(13)
      R1M5 = R1MACH(5)
      K = MIN(ABS(K1),ABS(K2))
      ELIM = 2.303E0*(K*R1M5-3.0E0)
      AB = AA*2.303E0
      ALIM = ELIM + MAX(-AB,-41.45E0)
      DIG = MIN(AA,18.0E0)
      FNUL = 10.0E0 + 6.0E0*(DIG-3.0E0)
      RL = 1.2E0*DIG + 3.0E0
      SLAK = 3.0E0+4.0E0*(-LOG10(TOL)-7.0E0)/11.0E0
      SLAK = MAX(SLAK,3.0E0)
      ERTOL = TOL*10.0E0**SLAK
      RM = 0.5E0*(ALIM + ELIM)
      RM = MIN(RM,200.0E0)
      RM = MAX(RM,RL+10.0E0)
      R2 = MIN(RM,FNUL)
      IF (KPRINT.GE.2) THEN
        WRITE (LUN,99998)
99998   FORMAT (' PARAMETERS'/
     *     5X,'TOL ',8X,'ELIM',8X,'ALIM',8X,'RL  ',8X,'FNUL',8X,'DIG')
        WRITE (LUN,99997) TOL, ELIM, ALIM, RL, FNUL, DIG
99997   FORMAT (1X,6E12.4/)
      ENDIF
C-----------------------------------------------------------------------
C     Set other constants needed in the tests.
C-----------------------------------------------------------------------
      CI = (0.0E0,1.0E0)
      HPI = 2.0E0*ATAN(1.0E0)
      RHPI = 1.0E0/HPI
      PI = HPI + HPI
C-----------------------------------------------------------------------
C     Generate angles for construction of complex Z to be used in tests.
C-----------------------------------------------------------------------
C     KDO(K), K = 1,IL  determines which of the IL angles in -PI to PI
C     are used to compute values of Z.
C       KDO(K) = 0  means that the index K will be used for one or two
C                   values of Z, depending on the choice of KEPS(K)
C              = 1  means that the index K and the corresponding angle
C                   will be skipped
C     KEPS(K), K = 1,IL determines which of the angles get incremented
C     up and down to put values of Z in regions where different
C     formulae are used.
C       KEPS(K)  = 0  means that the angle will be used without change
C                = 1  means that the angle will be incremented up and
C                   down by EPS
C     The angles to be used are stored in the T(I) array, I = 1,ITL.
C-----------------------------------------------------------------------
      IF (MQC.NE.2) THEN
        NL = 2
        IL = 5
        DO 5 I = 1,IL
          KEPS(I) = 0
          KDO(I) = 0
    5   CONTINUE
        KDO(5) = 1
        NUL = 5
        XNU(1) = 0.0E0
        XNU(2) = 1.0E0
        XNU(3) = 2.0E0
        XNU(4) = 0.5E0*FNUL
        XNU(5) = FNUL + 1.2E0
      ELSE
        NL = 4
        IL = 13
        DO 6 I = 1,IL
          KDO(I) = 0
          KEPS(I) = 0
    6   CONTINUE
        KDO(2) = 1
        KDO(6) = 1
        KDO(8) = 1
        KDO(11) = 1
        KDO(12) = 1
        KDO(13) = 1
        KEPS(3) = 1
        KEPS(4) = 1
        KEPS(5) = 1
        KEPS(9) = 1
        NUL = 6
        XNU(1) = 0.0E0
        XNU(2) = 0.6E0
        XNU(3) = 1.3E0
        XNU(4) = 2.0E0
        XNU(5) = 0.5E0*FNUL
        XNU(6) = FNUL + 1.2E0
      ENDIF
      I = 2
      EPS = 0.01E0
      FILM = IL - 1
      T(1) = -PI + EPS
      DO 30 K = 2,IL
        IF (KDO(K).EQ.0) THEN
          T(I) = PI*(-IL+2*K-1)/FILM
          IF (KEPS(K).NE.0) THEN
            TS = T(I)
            T(I) = TS - EPS
            I = I + 1
            T(I) = TS + EPS
          ELSE
            I = I + 1
          ENDIF
        ENDIF
   30 CONTINUE
      ITL = I - 1
C-----------------------------------------------------------------------
C     Test values of Z in -PI/2.lt.arg(Z).le.PI
C-----------------------------------------------------------------------
      IF (KPRINT.GE.2) THEN
        WRITE (LUN,99996)
99996   FORMAT (' CHECKS IN THE (Z,FNU) SPACE')
      ENDIF
      LFLG = 0
      DO 190 KODE = 1,2
        DO 180 N = 1,NL
          DO 170 NU = 1,NUL
C-----------------------------------------------------------------------
C     Construct values which will be used to set
C           COE1 = exp(i*(FNU+1)*PI/2) and
C           COE2 = (2/pi)*exp(-i*FNU*PI/2).
C-----------------------------------------------------------------------
            FNU = XNU(NU)
            IFNU = INT(FNU)
            FFNU = FNU - IFNU
            ARG = HPI*FFNU
            CSGN = CMPLX(COS(ARG),SIN(ARG))
            I4 = MOD(IFNU,4) + 1
            CSGN = CSGN*CIP(I4)
            CSPN = CONJG(CSGN)*CMPLX(RHPI,0.0E0)
            CSGN = CSGN*CI
            DO 160 ICASE = 1,3
              IRB = MIN(2,ICASE)
              DO 150 IR = IRB,4
C-------------- switch (icase)
                GO TO (50, 60, 70), ICASE
   50           CONTINUE
                  R = (EPS*(4-IR)+2.0E0*(IR-1))/3.0E0
                  GO TO 80
   60           CONTINUE
                  R = (2.0E0*(4-IR)+R2*(IR-1))/3.0E0
                  GO TO 80
   70           CONTINUE
                  IF (R2.GE.RM) GO TO 160
                  R = (R2*(4-IR)+RM*(IR-1))/3.0E0
   80           CONTINUE
C-------------- end switch
                DO 140 IT = 1,ITL
                  CT = COS(T(IT))
                  ST = SIN(T(IT))
                  IF (ABS(CT).LT.ATOL) CT = 0.0E0
                  IF (ABS(ST).LT.ATOL) ST = 0.0E0
                  Z = CMPLX(R*CT,R*ST)
                  XX = REAL(Z)
                  YY = AIMAG(Z)
                  CALL CBESI(Z, FNU, KODE, N, W, NZ2, IERR)
C---------------- Underflow in CBESI - skip test for this case.
                  IF (NZ2.NE.0) GO TO 140
                  CALL CBESK(Z, FNU, KODE, N, Y, NZ1, IERR)
C---------------- Underflow in CBESK - skip test for this case.
                  IF (NZ1.NE.0) GO TO 140
                  ZN = Z*CI
                  XN = -YY
                  YN = XX
                  CALL CBESY(ZN, FNU, KODE, N, V, NZ, CWRK, IERR)
C---------------- Underflow in CBESY - skip test for this case.
                  IF (NZ.NE.0) GO TO 140
                  COE1 = CSGN
                  COE2 = CSPN
                  IF (KODE.EQ.2) THEN
C------------------ Adjust scale for I and K functions.
                    CC = -XX - ABS(XX)
                    IF (CC.GT.(-ALIM)) THEN
                      CW = CMPLX(CC,-YY)
                      COE2 = COE2*EXP(CW)
                    ELSE
C-------------------- Scaling problem - skip test for this case
                      COE2 = CMPLX(0.0E0,0.0E0)
                      GO TO 140
                    ENDIF
                  ENDIF
                  DO 110 KK = 1,N
                    Y(KK) = Y(KK)*COE2
                    W(KK) = W(KK)*COE1
                    COE1 = COE1*CI
                    COE2 = -COE2*CI
  110             CONTINUE
C-----------------------------------------------------------------------
C     Compare Y(ZN,FNU) with COE1*I(Z,FNU)-COE2*K(Z,FNU).
C-----------------------------------------------------------------------
                  MFLG = 0
                  DO 120 I = 1,N
                    AB = FNU + I - 1
                    AA = MAX(0.5E0,AB)
                    CW = W(I) - Y(I)
                    AV = CABS(V(I))
                    ER = CABS(CW-V(I))
                    IF (AV.NE.0.0E0) THEN
                      IF (YN.EQ.0.0E0) THEN
                        IF (XN.GT.0.0E0) THEN
                          IF (ABS(XN).LT.AA) ER = ER/AV
                        ELSE
                          IF (ABS(FFNU-0.5E0).LT.0.125E0) THEN
                            IF (ABS(XN).LT.AA) ER = ER/AV
                          ELSE
                            ER = ER/AV
                          ENDIF
                        ENDIF
                      ELSE
                        ER = ER/AV
                      ENDIF
                    ENDIF
                    AER(I) = ER
                    IF (ER.GT.ERTOL) MFLG = 1
  120             CONTINUE
                  IF (MFLG.NE.0) THEN
                    IF (LFLG.EQ.0) THEN
                      IF (KPRINT.GE.2) THEN
                        WRITE (LUN,99995) ERTOL
99995                   FORMAT (/' CASES WHICH VIOLATE THE RELATIVE ',
     *                      'ERROR TEST WITH ERTOL = ',E12.4/)
                        WRITE (LUN,99994)
99994                   FORMAT (' INPUT TO CBESY   ZN, FNU, KODE, N')
                      ENDIF
                      IF (KPRINT.GE.3) THEN
                        WRITE (LUN,99993)
99993                   FORMAT (' COMPARE Y(ZN,FNU) WITH COE1*I(Z,FNU)',
     *                     '-COE2*K(Z,FNU)')
                        WRITE (LUN,99992)
99992                   FORMAT (' Z = ZN*EXP(-i*PI/2)'/
     *                     ' COE1 = EXP(i*(FNU+1)*PI/2)  ',
     *                     ' COE2 = (2/PI)*EXP(-i*FNU*PI/2)')
                        WRITE (LUN,99991)
99991                   FORMAT (' RESULTS FROM CBESY   V(KK)'/
     *                     9X,'FROM CBESI   W(KK)'/
     *                     9X,'FROM CBESK   Y(KK)')
                        WRITE (LUN,99990)
99990                   FORMAT (' TEST CASE INDICES   IR, IT, ICASE'/)
                      ENDIF
                    ENDIF
                    LFLG = LFLG + 1
                    IF (KPRINT.GE.2) THEN
                      WRITE (LUN,99989) ZN, FNU, KODE, N
99989                 FORMAT ('   INPUT:   ZN=',2E12.4,3X,'FNU=',E12.4,
     *                   3X,'KODE=',I3,3X,'N=',I3)
                    ENDIF
                    IF (KPRINT.GE.3) THEN
                      WRITE (LUN,99988) (AER(K),K=1,N)
99988                 FORMAT ('   ERROR:   AER(K)=',4E12.4)
                      WRITE (LUN,99987) Z, COE1, COE2
99987                 FORMAT (12X,'Z=',2E12.4/12X,'COE1=',2E12.4,3X,
     *                   'COE2=',2E12.4)
                      KK = MAX(NZ1,NZ2) + 1
                      KK = MIN(N,KK)
                      WRITE (LUN,99986) V(KK), W(KK), Y(KK)
99986                 FORMAT (' RESULTS:   V(KK)=',2E12.4/
     *                   12X,'W(KK)=',2E12.4/12X,'Y(KK)=',2E12.4)
                      WRITE (LUN,99985) IR, IT, ICASE
99985                 FORMAT ('    CASE:   IR=',I3,3X,'IT=',I3,3X,
     *                   'ICASE=',I3/)
                    ENDIF
                  ENDIF
  140           CONTINUE
  150         CONTINUE
  160       CONTINUE
  170     CONTINUE
  180   CONTINUE
  190 CONTINUE
      IF (KPRINT.GE.2) THEN
        IF (LFLG.EQ.0) THEN
          WRITE (LUN,99984)
99984     FORMAT (' QUICK CHECKS OK')
        ELSE
          WRITE (LUN,99983) LFLG
99983     FORMAT (' ***',I5,' FAILURE(S) FOR CBESY IN THE (Z,FNU) ',
     *       'PLANE')
        ENDIF
      ENDIF
      IPASS = 0
      IF (LFLG.EQ.0) THEN
        IPASS = 1
      ENDIF
      IF (IPASS.EQ.1.AND.KPRINT.GE.2) THEN
        WRITE (LUN,99982)
99982   FORMAT (/' ****** CBESY  PASSED ALL TESTS  ******'/)
      ENDIF
      IF (IPASS.EQ.0.AND.KPRINT.GE.1) THEN
        WRITE (LUN,99981)
99981   FORMAT (/' ****** CBESY  FAILED SOME TESTS ******'/)
      ENDIF
      RETURN
      END
